home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / DBio / DSeqFile.cpp < prev    next >
Text File  |  1996-07-05  |  12KB  |  488 lines

  1. // DSeqFile.cp 
  2. // d.g.gilbert
  3.     
  4.  
  5. #include <ncbi.h>
  6. #include <dgg.h>
  7. #include <stdio.h>
  8. #include <DFile.h>
  9. #include "DSeqFile.h"
  10. #include "DSeqList.h"
  11. #include "ureadseq.h"
  12. #include "DMegaSequence.h"
  13. #include "DSequence.h"
  14.  
  15.  
  16.  
  17. Global char* gOutputName = NULL;
  18.  
  19. Boolean DSeqFile::gWriteMasks = true;
  20. Boolean DSeqFile::gSaveMasks= true;
  21. static char* kMaskName = "#Mask";
  22.  
  23. // static
  24. void DSeqFile::DontSaveMasks()
  25. {
  26.     gSaveMasks= gWriteMasks;
  27.     gWriteMasks= false;
  28. }
  29.  
  30. void DSeqFile::DoSaveMasks()
  31. {
  32.     gWriteMasks= gSaveMasks;
  33. }
  34.  
  35.  
  36. Boolean dupBases = false; // !?
  37.  
  38. //static
  39. DSequence* DSeqFile::MakeSequence(char* name, char*& bases, char* info, long modtime)
  40. {
  41.     char *b;
  42.     long nbases;
  43.     Boolean ismega;
  44.     DSequence* aSeq;
  45.     
  46.     for (b= bases, nbases=0; 
  47.             (b) && (*b) && nbases<DMegaSequence::kMegaPartSize; 
  48.             b++, nbases++)
  49.                 ;
  50.     ismega= nbases>=DMegaSequence::kMegaPartSize; 
  51.      
  52.     if (ismega) {
  53.         aSeq= new DMegaSequence();
  54.         } 
  55.     else {
  56.         aSeq= new DSequence();
  57.         }
  58.     aSeq->SetName(name);
  59.     aSeq->SetBases(bases, dupBases);
  60.     aSeq->SetInfo(info);
  61.      if (modtime!=0) aSeq->SetTime( (unsigned long)modtime);
  62.     aSeq->UpdateFlds();      
  63.     return aSeq;
  64. }
  65.  
  66.  
  67.  
  68. #if NEED_LATER
  69. // must be in scope of a DSeqFile:: method
  70. Local const  struct formatTable {
  71.   char  *name;
  72.   short num;
  73.   } formname[] = {
  74.     {"ig",  kIG},
  75.     {"stanford", kIG},
  76.     {"genbank", kGenBank},
  77.     {"gb", kGenBank},
  78.     {"nbrf", kNBRF},
  79.     {"embl", kEMBL},
  80.     {"gcg", kGCG},
  81.     {"uwgcg", kGCG},
  82.     {"dnastrider", kStrider},
  83.     {"strider", kStrider},
  84.     {"fitch", kFitch},
  85.     {"pearson", kPearson},
  86.     {"fasta", kPearson},
  87.     {"zuker", kZuker},
  88.     {"olsen", kOlsen},
  89.     {"phylip", kPhylip},
  90.     {"phylip3.2", kPhylip2},
  91.     {"phylip3.3", kPhylip3},
  92.     {"phylip3.4", kPhylip4},
  93.     {"phylip-interleaved", kPhylip4},
  94.     {"phylip-sequential", kPhylip2},
  95.     {"plain", kPlain},
  96.     {"raw", kPlain},
  97.     {"pir", kPIR},
  98.     {"codata", kPIR},
  99.     {"asn.1", kASN1},
  100.     {"msf", kMSF},
  101.     {"paup", kPAUP},
  102.     {"nexus", kPAUP},
  103.     {"pretty", kPretty},
  104.     {0,0}
  105.   };
  106. #endif
  107.  
  108. enum { kMaxFormatNames = 18 }; // != kMaxFormat ??
  109. static char *gSeqFileFormats[] = { //DSeqFile::kMaxFormat+1
  110.         "IG|Stanford",
  111.         "GenBank|GB",
  112.         "NBRF",
  113.         "EMBL",
  114.         "GCG",
  115.         "DNAStrider",
  116.         "Fitch",
  117.         "Pearson|Fasta",
  118.         "Zuker [in only]",
  119.         "Olsen [in only]",
  120.         "Phylip3.2|Phylip2",
  121.         "Phylip|Phylip4",
  122.         "Plain|Raw",
  123.         "PIR|CODATA",
  124.         "MSF",
  125.         "ASN.1",
  126.         "PAUP|NEXUS",
  127.         "Pretty [out only]",// # 18
  128.         //"AutoSeq [in only]", 
  129.         0};
  130.  
  131. const char* DSeqFile::FormatName(short index)
  132.     if (index<1 || index>kMaxFormatNames) return "";
  133.     else return gSeqFileFormats[index-1];
  134. }
  135.  
  136. short  DSeqFile::FormatFromName(const char* name)
  137. {
  138.     if (!name) return kUnknown;
  139.     long   i, len1, namelen= StrLen(name);
  140.     char      *form1, *form2;
  141.     if (namelen<2) return kUnknown;
  142.     for (i=0; i<kMaxFormat; i++) {
  143.         form1= (char*) gSeqFileFormats[i];
  144.         form2= StrChr(form1,'|'); 
  145.         if (form2) { len1= form2 - form1; form2++; }
  146.         else len1= StrLen( form1);
  147.         if (StrNICmp(name, form1,len1) == 0) return i+1;
  148.         else if (form2 && StrICmp(name, form2) == 0) return i+1;
  149.         }
  150.     return kUnknown;
  151. }
  152.  
  153.  
  154. char* DSeqFile::FixSeqID( char* s, short leftmargin)
  155. {
  156.     register char * cp, * ep;
  157.     for (cp= s; *cp && isspace(*cp); cp++) ;
  158.     for (ep= cp; *ep && !isspace(*ep); ep++) ;
  159.     if (ep>cp && ep[-1] == ',') ep--; // fix for fasta "name, #bases, # checksum"
  160.     if (*ep) *ep= 0;
  161.     return cp;
  162. }
  163.  
  164. short DSeqFile::SeqFileFormatWrapper( DFile* aFile) 
  165. {
  166.     long    skiplines = 0;
  167.     short    aFormat, err = 0;
  168.  
  169.     aFile->Open("r");  
  170.     aFormat= seqFileFormatFp( aFile->fFile, &skiplines, &err);
  171.     aFile->Close();  
  172.     return    aFormat;
  173. }
  174.  
  175.  
  176. short DSeqFile::ReadSeqFile( DFile* aFile, DSeqList* aSeqList) 
  177. {
  178.     char    seqid[256], *seqidptr, *bases, *seqlist;
  179.     long    skiplines = 0, skipi, seqlen;
  180.     short aFormat, atseq = 0, whichSeq = 0, nseq = 0, err = 0;
  181.     char  *defname;
  182.     Boolean notdone = true, needrewind = false;
  183.     DSequence* aSeq;
  184.     
  185.     aFile->OpenFile("r");    
  186.     aFormat = seqFileFormatFp( aFile->fFile, &skiplines, &err);
  187.     aFile->OpenFile(); // rewind
  188.     defname= FixSeqID( gDefSeqName, 0);
  189.     
  190.     switch (aFormat) {
  191.         case kUnknown:
  192.         case kNoformat:
  193.             aFile->CloseFile();
  194.             return aFormat;
  195.          
  196.             // all multi-pass (interleaved) formats go here !        
  197.         case kOlsen:
  198.         case kMSF:
  199.         case kPAUP:
  200.         case kPhylip2:
  201.         case kPhylip4:
  202.             needrewind = true;
  203.             aFile->CloseFile();
  204.             seqlist = listSeqs( aFile->GetName(), skiplines, aFormat, &nseq, &err);
  205.             MemFree(seqlist); 
  206.             break;
  207.         }
  208.         
  209.     skipi= skiplines;
  210.     while (notdone) {
  211.             // this one-pass thru gFileRef may well be buggy for current readSeq
  212.             // need to rewind file & count whichSeq++ for each sequence...?
  213.         if (needrewind) {
  214.             whichSeq++; 
  215.             atseq= 0; 
  216.             aFile->OpenFile("r");
  217.             skipi= skiplines;
  218.             }
  219.         else {
  220.             whichSeq = 1; 
  221.             atseq= 0;
  222.             }
  223.         seqlen = 0; seqidptr= seqid; *seqidptr= 0;
  224.         bases = readSeqFp( whichSeq, aFile->fFile, skipi, aFormat, &seqlen, 
  225.                                                         &atseq, &err, seqidptr);
  226.         skipi= 0; 
  227.         
  228.                 // install bases into new TSequence, and stick TSequence into list
  229.         if (bases) {
  230.             if (seqlen > 0) {
  231.                 char *name, * cp;
  232.                 Boolean ismask= false;
  233.                 
  234.                 if (*seqidptr != 0) {
  235.                     if ( (cp= StrStr(seqidptr, kMaskName)) != NULL) {
  236.                         ismask= true;
  237.                         *cp= 0;
  238.                         }
  239.                     name= FixSeqID( seqidptr, 0);
  240.                     }
  241.                 else 
  242.                     name= defname;
  243.                     
  244.                 if (ismask) {
  245.                     aSeq= aSeqList->FindNamedSeq( name);
  246.                     if (aSeq) {
  247.                         // convert mask...
  248.                         register char *mp;
  249.                         for (mp= bases; *mp; mp++) { 
  250.                             *mp -= '?'; // to zero relative
  251.                             *mp |= 0x80; // add high bit
  252.                             }
  253.                         aSeq->SetMasks( bases);
  254.                         }
  255.                     }
  256.                 else {
  257.                     aSeq = MakeSequence( name, bases, NULL, 0); //gFileRef->moddate
  258.                     aSeqList->InsertLast(aSeq);
  259.                     }
  260.                 }
  261.             else notdone= false;
  262.             MemFree(bases); // readSeq made this ptr w/ malloc call... MakeSeq copies this data... 
  263.             }
  264.         else notdone= false; // nodata -- if format=0, no data is read
  265.             
  266.         if (needrewind) notdone &= (whichSeq < nseq); // stop now
  267.         else notdone &= !aFile->Eof();  
  268.         }
  269.         
  270.     aFile->CloseFile();
  271.     return aFormat;
  272. }
  273.  
  274.  
  275.  
  276.  
  277. void    DSeqFile::WriteSeqHeader( DFile* aFile,  
  278.                                 short format, short nseqs, long nbases, char* outputName,
  279.                                 Boolean& needSameSize, Boolean& isInterleaved)
  280. {
  281.     const char *kASN1headline = "Bioseq-set ::= {\nseq-set {\n";
  282.     char    line[256], name[64];
  283.     
  284.     // this is a ureadseq.h call
  285.     // need to do somewhere once, then let user set prefs...
  286.  
  287. #define gPrettyInit1(p) { \
  288.   p.isactive=false;\
  289.   p.baseonlynum=true;\
  290.   p.numline= p.atseq= 0;\
  291.   p.numright= p.numtop= true; \
  292.   p.numleft= p.numbot= false;\
  293.   p.nameright= true; \
  294.   p.nameleft= p.nametop= false;\
  295.   p.noleaves= p.domatch= p.degap= false;\
  296.   p.matchchar='.';\
  297.   p.gapchar='-';\
  298.   p.namewidth=8;\
  299.   p.numwidth=5;\
  300.   p.interline=1;\
  301.   p.spacer=10;\
  302.   p.seqwidth=50;\
  303.   p.tab=0; }
  304.  
  305.     gPrettyInit1(gPretty);
  306.     //gPretty.interline= 2; // not active, along w/ numtop,numbot,nametop...
  307.     // need to work into this class
  308.     
  309.     
  310.         
  311.     needSameSize= false;
  312.     isInterleaved= false;
  313.     aFile->Open("a"); 
  314.  
  315.     if (outputName) StrNCpy( name, outputName, 63);
  316.     else StrCpy( name, DFileManager::kUntitled);
  317.     
  318.     switch (format) {
  319.         case kPIR:
  320.             aFile->WriteLine( "\\\\\\\n");
  321.             break;
  322.         case kASN1:
  323.             aFile->WriteLine( (char*)kASN1headline);
  324.             break;
  325.  
  326.         case kMSF:     // this is for top of interleaved file
  327.             {
  328.             isInterleaved= true;
  329.          unsigned long checksum = 0, checkall = 0;
  330.             //checksum= GCGchecksum(seq, seqlen, &checkall);
  331.             sprintf(line,"\n %s  MSF: %d  Type: N  June 01, 1776 12:00  Check: %d ..\n\n",
  332.                                         name, nbases, checkall);
  333.             aFile->WriteLine( line);
  334.             }
  335.             break;
  336.  
  337.         case kPhylip2:     
  338.             isInterleaved= false;
  339.             goto allphyl;
  340.             break;
  341.  
  342.         case kPhylip4:      
  343.             isInterleaved= true;
  344.         allphyl:
  345.             {
  346. //    if (phylvers >= 4) fprintf(foo," %d %d\n", seqout, seqlen);
  347. //    else fprintf(foo," %d %d YF\n", seqout, seqlen);
  348.             if (isInterleaved) sprintf( line, " %d %d\n", nseqs, nbases);  // " %d %d F\n"
  349.             else sprintf( line, " %d %d I \n", nseqs, nbases);  // " %d %d FI \n"
  350.             aFile->WriteLine( line);
  351.             needSameSize= true;
  352.             }
  353.             break;
  354.  
  355.         case kPAUP:
  356.             isInterleaved= true;
  357.             needSameSize= true;
  358.             aFile->WriteLine( "#NEXUS\n");
  359.             sprintf( line,"[%s -- data title]\n\n", name);
  360.             aFile->WriteLine( line);
  361.             break;
  362.  
  363.         case kPretty:
  364.             isInterleaved= true;
  365.             gPretty.isactive= true;
  366.             break;
  367.             
  368.         default:
  369.             break;
  370.         }
  371.         
  372.     aFile->CloseFile();
  373. }
  374.  
  375.  
  376.  
  377. void    DSeqFile::DoInterleave( DFile* toFile, DFile* fromFile, 
  378.                                 short outform, short seqtype, short nseqs, long nbases,
  379.                                 short nlines, short noutindex, long* outindex)
  380. {
  381.     short        leaf, iline, iseq;
  382.     char        *cp, line[256];
  383.     
  384.     toFile->Open("a");
  385.     fromFile->Open("r");
  386.     
  387.     for (leaf=0; leaf<nlines; leaf++) {
  388.     
  389.         if (outform == kMSF && leaf == 1) { // do after writing seq names
  390.             toFile->WriteLine( "//\n\n");
  391.             }
  392.             
  393.         if (outform == kPAUP && leaf==1) { // do after writing seq names
  394.             switch (seqtype) {
  395.                 case DSequence::kDNA     : cp= "dna"; break;
  396.                 case DSequence::kRNA     : cp= "rna"; break;
  397.                 case DSequence::kNucleic : cp= "dna"; break;
  398.                 case DSequence::kAmino   : cp= "protein"; break;
  399.                 case DSequence::kOtherSeq: cp= "dna"; break;
  400.                 }
  401.             toFile->WriteLine( "\nbegin data;\n");
  402.             sprintf(line," dimensions ntax=%d nchar=%d;\n", nseqs, nbases); 
  403.             toFile->WriteLine( line);
  404.             sprintf(line," format datatype=%s interleave missing=%c", cp, gPretty.gapchar); // gapchar == indelHard ??? 
  405.             toFile->WriteLine( line);
  406.             //if (gPretty.domatch) { sprintf(line," matchchar=%c", gPretty.matchchar); 
  407.                 //Mwrite( toFile, line); }
  408.             toFile->WriteLine( ";\n  matrix\n");
  409.             }
  410.  
  411.         for (iseq=0; iseq<noutindex; iseq++) {
  412.             fromFile->SetDataMark( outindex[iseq]);
  413.             
  414.             for (iline=0; iline<=leaf; iline++)
  415.                 if (fromFile->ReadLine( line, 256) != 0) *line= 0;
  416.                 
  417.             if (fromFile->Tell() <= outindex[iseq+1])
  418.                 toFile->WriteLine( line, false);  //remember- ReadLine a.k.a fgets retains newline
  419.             }
  420.  
  421.      for (iline=0; iline<gPretty.interline; iline++)
  422.         toFile->WriteLine("\n"); //LineEnd 
  423.             
  424.         }
  425.         
  426.     fromFile->Close();
  427.     toFile->Close();
  428. }
  429.  
  430.  
  431. void    DSeqFile::WriteSeqTrailer( DFile* aFile,
  432.                                 short format, short /*nseqs*/, short /*nlinesperseq*/, long /*nbases*/)
  433. {
  434.     switch (format) {
  435.   
  436.       case kASN1:  
  437.             aFile->Open("a"); 
  438.             aFile->WriteLine("} }\n"); 
  439.             aFile->Close(); 
  440.             break;
  441.     
  442.       case kPAUP:    
  443.             aFile->Open("a"); 
  444.             aFile->WriteLine(";\n  end;\n");
  445.             aFile->Close(); 
  446.             break;    
  447.         }
  448. }
  449.  
  450.  
  451.  
  452. void DSeqFile::WriteSeqWrapper( DFile* aFile,
  453.                                     char *seq, long seqlen, short outform, char *seqid)
  454.     if (!aFile->IsOpen())  aFile->Open("a"); 
  455.     seqid= StrDup(seqid);
  456.     char* cp= FixSeqID(seqid, 0);
  457.     gLinesPerSeqWritten= writeSeq( aFile->fFile, seq, seqlen, outform, cp);
  458.     MemFree(seqid);
  459.     //aFile->Close(); //?? want left open for appending by caller ? e.g., seqmail
  460. }  
  461.  
  462. void DSeqFile::WriteMaskWrapper( DFile* aFile,
  463.                                     char *mask, long seqlen, short outform, char *seqid)
  464.     char maskid[128];
  465.     
  466.     if (!aFile->IsOpen())  aFile->Open("a"); 
  467.     seqid= StrDup(seqid);
  468.     char* cp= FixSeqID(seqid, 0);
  469.     StrNCpy( maskid, cp, sizeof(maskid)-5);
  470.     StrNCat( maskid, kMaskName, sizeof(maskid));
  471.     MemFree(seqid);
  472.         //     massage mask into writable form
  473.     mask= StrDup(mask); 
  474.     Boolean hasdata= false;
  475.     for (cp= mask; *cp; cp++) { 
  476.         *cp &= 0x7F; // drop hi bit
  477.         if (*cp != 0) hasdata= true;
  478.         *cp += '?';  // set 0 == '?', mask range is mostly 'A-Za-z'
  479.      }
  480.     if (hasdata)
  481.         gLinesPerSeqWritten= writeSeq( aFile->fFile, mask, seqlen, outform, maskid);
  482.     MemFree(mask);
  483. }  
  484.  
  485.